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We compute the thermodynamic properties of the glass phase in a bi- 
nary mixture of soft spheres. Our approach is a generalization to mixtures 
of the replica strategy, recently proposed by Mezard and Parisi, providing 
a first principle statistical mechanics computation of the thermodynamics of 
glasses. The method starts from the inter-atomic potentials, and translates 
the problem into the study of a molecular liquid. We compare our analytical 
predictions to numerical simulations, focusing onto the values of the thermo- 
dynamic transition temperature and the configurational entropy. 
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I. INTRODUCTION 



In this paper we present the generahzation to the binary mixture case of a thermodynamic 
theory of glasses, recently proposed in which allows to deduce equilibrium properties 
of fragile glasses [Q] from those of the corresponding liquid phase, computed for a molecular 
liquid consisting of m 'clones' of the system with m < 1. 

The hypothesis at the heart of this strategy is the existence of a liquid-glass thermody- 
namic transition, driven by the 'entropy crisis' predicted by Kauzmann P], and the scenario 
is similar to the one described by Adam, Gibbs and Di Marzio P-Q. The transition con- 
sidered here can be also explained in terms of a certain type of replica symmetry breaking 
(called 'one step replica symmetry breaking' -IRSB). It shares its main features with the 
glass transition found in some discontinuous spin-glasses model, as first proposed by Kirk- 
patrick Thirumalai and Wolynes ||^. 

We identify the mode coupling temperature Tmct Q with the dynamical temperature 



Td which exists in discontinuous spin-glasses [|rT]-[T^, and we assume that below this tem- 
perature the phase space can be partitioned in a very large number of different free energy 
valleys. These valleys are supposed to be, in terms of free-energy, the equivalent of the 
so called inherent structures [|14[, which are built from the minima of the potential energy 
together with their basins of attraction. 

In other words, we suppose that, for T < Tmct, a typical equilibrium configuration 
belongs to one of these valleys. We label the valleys with an index a, and denote for each 
valley the free energy density as /„, the subset of equilibrium configurations belonging to 
the valley as Vq, and the corresponding restricted partition function as Z^. The canonical 
partition function can then be written in the following way: 

Z ^Y.Za = I dxe-"""^^^ = Y^e-^'^f- (1) 

a a 

where the function B. is the Hamiltonian of the system and (5 is the inverse temperature. The 
number of valleys with a given value of free energy density is defined as A/'(/) = ^{f~fa), 
and we assume that in the thermodynamic limit it becomes a continuous function. It is then 
possible to write the partition function as: 

df Ar(/)e-^^^ = fdf e-^[/3/-s(/,T)] 



where we have introduced the complexity S = \ogAf/N. 

Let us note that the system in equilibrium does not minimize the free energy of the single 
valleys, but a 'collective' thermodynamic potential 0(T), that we interpret as the actual free 
energy in the liquid and glass phases. 0(T) is defined by: 

0(r) = /*-TS(r,r). (3) 

where /* is the temperature dependent free energy which minimizes the function / — 

TS(/,r). 

In this picture, the total entropy density is the sum of the entropy inside the valley, and 
of the entropy coming from the very large number of valleys, with the same value of free 
energy, that the system is allowed to explore: 
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Stot — S valley + S(T, /*) , 



(4) 



hence the complexity introduced here is completely equivalent to the usual concept of con- 
figurational entropy of a super-cooled liquid. 

Assuming the existence of this decomposition of phase space into valleys, we will show 
that there is a finite temperature (the so-called Kauzmann temperature) where the 
system undergoes a thermodynamic transition with the following features: 

• Tk is the temperature where the complexity S vanishes l]!],^. This means that, unlike 
the liquid phase, in the whole low temperature glass phase, only a non exponential 
number of valleys contribute to the partition function, namely the ones with the lowest 
free energy density fmin- 

• At Tk there is a second order transition from the thermodynamic point of view. The 
free energy is continuous and there is no latent heat. The specific heat 'jumps' from 
the liquid value to a smaller one, in agreement with the Dulong and Petit law. 

• At Tk there is a discontinuity of the order parameter. Below Tk-, in the glass phase, 
the system is an amorphous solid and the thermal average of the local particle density 
becomes non uniform, exhibiting peaks at the favoured positions where the particles 
tend to be trapped in some cages. The order parameter is related to the spatial mod- 
ulation of the density, and it goes discontinuously from zero in the high temperature 
liquid phase to a finite modulation in the glass phase. 

This transition could be experimentally observed only if one would be able to cool the 
liquid at an infinitely slow rate, and should correspond to the temperature where the 
viscosity is supposed to diverge (following for instance a generalized Vogel-Fulcher law 7] oc 
exp(T — Tk)^^) 1^. In real experiments, infinitely slow cooling is not available, and the 
correlation time becomes of order of the experimental time at a temperature T^, which is in 
general an intermediate temperature Tk <Tg < Tmct- The value of Tg could be computed 
only if we had under control the time dependence of the correlation functions. In this paper 
we study only static quantities, and we cannot say anything about the value of Tg or the 
temperature dependence of the viscosity above Tk- 

In such an 'entropy crisis' scenario, it has been shown [p!5|,3 that the thermodynamic 
properties of the glass phase can be computed in principle by considering m replicas of the 
original system, constrained to stay in the same valley, by means of a small but extensive 
coupling term. In this case, the arguments used in the derivation of (|^) can be applied again, 
leading to a replicated version of the same equation: 

$(m,T) =Min\f {m f - TJ:{f,T)) . (5) 

Once again, each of the m systems does not reach the lowest possible free energy, but the 
one which optimizes the balance between the free energy and complexity in (^). 

Interestingly enough, one can derive many properties of the system from (^ if one is 
able to continue it analytically and compute it for any real value of m, thinking about m 
as a new parameter of the problem. Indeed, $(m, T), considered as a function of m, gives 
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access to the configurational entropy (complexity) T) through a Legendre transform. 
This imphes the relations: 

^_(9$(m,T) ^_m'^dmm,T)/m) 
dm T dm 

from which it is possible to eliminate m, obtaining which measure the number of 

valleys with a given value of free energy /. Let us underline that (||) gives access to the 
full curve of complexity versus free energy, while the equilibrium free energy of the physical 
system is obtained only after taking the limit m — > 1. 

As we shall see, the thermodynamic potential $(m)/m is a convex function of m with a 
maximum at a point m*{T), which is an increasing function of T, vanishing at T = 0. The 
second equation of (^) is thus well-defined for m < m*{T). At m = m*{T), the resulting 
complexity vanishes S = and the free energy d^/dm reaches a value fmm- For m < m*{T) 
the complexity S is non zero: it is thermodynamically favourable to select some valleys which 
have a free energy density larger than fmin, because of the corresponding gain in complexity. 
If one increases m beyond m*{T), the formula gives an unphysical negative complexity. 
In fact in the whole region m > m*{T) the correct value of / is / = fmin, and the complexity 
is zero. 

This is easily understood from the physical interpretation of the transition which we now 
turn to. The above scenario has: 

• a high temperature phase where m*(T) > 1. In this phase, when the limit m — 1 is 
performed and the equilibrium free energy of the original system is recovered, one gets 
a value feq > fmin together with a positive configurational entropy. 

• a low temperature phase where m*{T) < 1. In this phase, in the limit m — > 1 the 
equilibrium free energy is fmin and the configurational entropy is null. 

It is quite easy, at this point, to recognize these two thermodynamic phases as the super- 
cooled liquid one (high) and the glass one (low), separated by a thermodynamic transition of 
second order, driven by the vanishing of the configurational entropy, at the temperature 
where m*{TK) = 1- All the thermodynamic quantities in the glass phase can be computed 
from the replicated free energy (|^) at the point m*, which play the role of the free energy 
of the glass. 

This scenario of the glass transition is identical to the phase transition appearing in 
discontinuous (IRSB) spin glasses where it was first explained [|16[. The simplest example of 



such a discontinuous spin glass transition is the Random Energy Model [T^ which displays 
a total freezing at T = Tk- As first noticed by Kirkpatrick Thirumalai and Wolynes 0, 
discontinuous spin glasses provide some well defined mean field systems where the old ideas 
of Adam-Gibbs-Di Marzio of a real thermodynamic transition driven by entropic reasons are 
at work. The present approach allows to apply the replica method directly to the structural 
glasses (in spite of the absence of any quenched disorder in the Hamiltonian) . Assuming 
that the structural glass transition is characterized in the replica language by a IRSB, as 



in discontinuous spin glasses p8|-p0|, we can compute the thermodynamic properties of the 
glass phase. The comparison with the numerical results allows then to justify a posteriori 
the main hypothesis. 
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At this stage, m appears as an auxiliary parameter which may be interpreted as the 
effecive temperature of the vallyes; moreover one finds that in the low temperature phase 
1 — m gives the probability of finding two systems in the same valley [|1^ . 

Summarizing, the study of the liquid-glass transition and the investigation of the low 
temperature phase can be accomplished by computing the free energy of a replicated system 
in its liquid phase, or in other words 0, the free energy of a molecular liquid where each 
molecule has m atoms. The thermodynamic properties of the glass phase can be deduced 
by means of the analytic continuation to arbitrary real values of this parameter. 

In the previous works |]T],|^ this general approach was applied to a pure soft sphere sys- 
tem. The extension to binary mixtures is particularly important since there are well known 
examples of glass forming binary mixtures where an appropriate choice of the interaction 
parameters strongly inhibits crystalization. This allows therefore to get numerical results 
which can be compared to the analytical ones. Here we will consider in particular a mixture 
of soft spheres. 

After discussing the model in sect. II, we will present in sect. Ill the generalization to 
binary mixtures both of the small cage expansion and of the harmonic re-summation scheme 
introduced previously to deal with the molecular fiuid. Sect. IV describes the application 
of the HNC approximation to the center of mass degrees of freedom of the molecular fiuid. 
In the last section we will discuss our analytic results, together with some strategies for 
evaluating numerically the glass transition temperature Tr- and the configurational entropy 
behavior, and a comparison between the numerical estimates and those obtained analytically. 



II. GENERAL FRAMEWORK 

We study mixtures composed of two types of particles called + and — , with pairwise 
interactions. The Hamiltonian of our problem is: 

ff= ^ V^^^.^a:,-x,) e, €{-,+}, (7) 

l<i<j<Ar 

where the N particles move in a volume of a (i-dimensional space, and V'^'^, V'^^ , V 
are arbitrary short range interaction potentials. We call c+ (resp. c_) the fraction of -|- 
(resp. — ) particles. 

In the explicit computations described in the next section, we have chosen a binary 
mixture of soft spheres that has been extensively studied in the past through numerical 
simulations [ P^pT| . The potentials are given by: 

V''(r)=(^y\ (8) 



where 



1.2, a+_ = . (9) 



a 



The concentration is taken as c+ = 1/2, and the choice of the ratio R = a^+/a = 1.2 

is known to strongly inhibit crystalization. We also make the usual choice of considering 
particles with average diameter 1 by setting 
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(a++)3 + 2(a+_)3 + (a__)3 ^ ^ 
4 

All thermodynamic quantities depend on the density p = N/V and temperature T only 
through the parameter T = pT~^/^. For F larger than F^i = 1.45 (corresponding to lower 
temperatures) the dynamics becomes very slow and the autocorrelation time is very large. 
Hence the system enters the 'aging' regime, where violations of the equilibrium fluctuation- 



dissipation theorem are observed [20|. This value of F^) is supposed to correspond to the 



mode coupling transition below which the relaxation is dominated by activated processes 



21| . If this simple model behaves like a real fragile glass the Kauzmann transition, charac- 
terized by a discontinuity in the specific heat, is located below the dynamical transition, and 
cannot be directly accessed by numerical simulations, maybe with the exception of studies 



done on very small samples . 



The application of the theory to a more realistic potential, namely a Lennard- Jones 
binary mixture, will be treated in detail in a forthcoming paper . 



As previously explained, in order to obtain some information about the super-cooled 
liquid-glass thermodynamic transition, we consider the thermodynamics of a molecular liq- 
uid, whose molecules are composed of m atoms, each carrying a different replica index. 
The tendency to form molecules is forced by a small but extensive coupling term between 
particles of different replicas @]. Unlike the pure case, we are dealing here with a situation 
where particles are not all indistinguishable: we have particles of the " +" type and of the 
"— " type. Physically this has an important effect when R is not close to one. At i? ~ 1, 
it is clear that the valleys of the mixture are close to those of the pure system. More pre- 
cisely, taking one given valley of the pure {R = 1) system, one can generate A^!/A^+! A^_! 
valleys of the mixture with i? ~ 1, by choosing at random the positions of the + and the 
— particles: in this limit the main effect of the mixture is to add a factor to the entropy, 
whose value is A^log2 when N± = N/2. On the other hand, when R is very different from 
one, the valleys of the mixture are very different from those of the pure system; one cannot 
find a new valley by just exchanging a + particle with a — particle. This physical problem 
has an exact counterpart in replica space. One could study the case where molecules are 
formed by one particle of each of the m different replicas, irrespective of their ± nature. 
Qualitatively speaking, this would mean that interchanging two particles of different types, 
the two replicas to which particles belong would remain in the same valley, that is their 
free-energy would not change. There are two extreme possibilities, corresponding to the two 
cases discussed above: 



• For R very near to one, the system behaves similarly to the system at i? = 1. One can 
form molecules with particles of any type, and the exchange of a "+" particle with a 
"— " one gives a very small change in free energy. 

• For R quite different from 1, the exchange of a "+" particle with a "— " one is a process 
that can be safely neglected, since it gives a variation in energy that is much larger 
than kT. In this second case the molecules are built up of atoms of the same type. 

In each of these extreme cases the computation is simple: in the first case it just reduces 
to the R = 1 computation. In the second case, we can assume, as we shall do here, that each 
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molecule is built from m atoms which are all of the same type (all " +" or all " — " ). Then one 
only needs considering attractive coupling terms only between particles of the same kind. 
The computations in the crossover region are rather complex. For our case R = 1.2, we have 
decided to neglect this kind of corrections and to consider the molecules consisting only of 
particles of the same type. 

The replicated partition function is: 

- ^,4.^ EE / n^'^-" exp (-^ E - 4)+ 

JV + . JV_. J y z, ^^.^^ 

- E E^+(^<^4i) -^<^6«) - E E^-(^T„(i)-a;n(i)) h (11) 

ie{+} 0,^6 ie{-} 0-+'^ J 

where the sum over permutations of atoms in each molecule is taken into account, and 
= c+A^, = c^N . When relabeling particles, so that particle i of a given type in 
replica a corresponds to particle i of the same type in replica h (which is supposed to belong 
to the same molecules) and so on, the sum over permutations gives a factor (iV+! A^_!)*^"^~-^). 

As discussed in the previous section, in the glass phase the replicas becomes correlated, 
so the study of the transition is accomplished by choosing as order parameters the m-points 
correlation functions for each of the two different types of particles: 

p+(r\...,r-)= <?>{x\-r^)...b{xr-r"^)> . (12) 



p_(r\...,0= <?>{A-r^)...K^r-n> ■ (13) 
ie{-} 

The transition is signaled, then, by the onset of an off-diagonal non trivial correlation in 
replica space at T^, when the coupling functions ujj^ are sent to zero. This feature is studied 
as usual introducing the Legendre transform of the molecular (replicated) free energy: 

G\p\ = hm - logZ^[a;]- - / Hd'r'^ Y P.(r\ r™)iy,(r\ ^ (14) 

777, —> 1 ^ o=l £=+,- 

N ^ oo 



with 

W,{r\...,r'^) = Y^e{r'' -r^) (15) 

a<b 

Performing the limit u± — > is equivalent to searching a saddle point of the functional G[p]. 
In the presence of a glassy transition we expect the following behavior of order parameters 
and thermodynamic quantities: 

• For T > Tk the free energy is the liquid one (m = 1) and the order parameters are 
trivial, i.e. p±(r^) = c±p. 
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• For T < Tk, the glass free energy is the maximum with respect to m of the replicated 
free energy, and this maximum is found at m* < 1. The correlations p± become non 
trivial. From the free energy at the maximum we can compute all the thermodynamic 
quantities. 

The free energy and his first derivatives are continuous at T^, while the heat capacity falls 
suddenly from liquid-like to solid-like values when the temperature is decreased through T^- 
The transition, then, is of second order from the point of view of thermodynamics, but it 
is discontinuous in the order parameter which abruptly becomes a non trivial function of 
positions in different replicas. 

It is natural to describe the particle positions in term of center of mass coordinates 
and relative displacements m° with x1 = Zi + u°; and J2a u°; = 0. A useful simplification is 
the choice, for the polarising potentials uj±, of a quadratic coupling that allows to rewrite 
m) as: 



\i=l / \a=l 1=1 



N 



nh'^(E<) 



.i=l 



a=l 



exp -/3 E E ^"^^^ - h «I) + 

\ a=l i<j 



(16) 



In the absence of the interacting potential V ^ the {u\^ for a given % are Gaussian random 
variables with a vanishing first moment and a second moment given by 



1 

m 



(17) 



III. REPLICATED FREE-ENERGY 

A. Harmonic re-summation 

We are interested in the regime of low temperatures, where the molecules are expected to 
have a small radius, justifying a quadratic expansion of V in the partition function (|16|). 
After integrating over these quadratic fluctuations, one obtains: 

= J n d\ exp -Pm Y: V-^^^ {z, -z,)- "^Tt log {(3M) (18) 

1=1 \ i<j J 

where the matrix M, of dimension Nd x Nd, is given by: 

^K.^) = (e - + r^) - (^^ - ^) (19) 
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and Vfj^yir) = d'^v/dr^dvu (the indices /i and z/, running from 1 to d, denote space directions). 
We have thus found an effective Hamiltonian for the centers of masses Zi of the molecules, 
which basically looks like the original problem at the effective temperature T* = l/{(3m), 
complicated by the contribution of vibration modes. We shall proceed by using the same 
set of approximations which was proposed in the previous papers |1|,0]. We first perform 
a 'quenched approximation', which amounts to neglecting the feedback of vibration modes 
onto the centers of masses, substituting thus the term Tr log (/5M) in (^) by its expectation 
value, for center of mass positions Zi equilibrated at the temperature T*. This approximation 
becomes exact close to the Kauzmann temperature where m — > 1. 

Let us introduce the mean values of the diagonal terms of the matrix M: 

r, = 5: CP / d'rg^^'{r)]-AV^'' + - , (20) 

where the g^'^'{r) are the pair correlation functions. We neglect the fluctuation of these 
diagonal terms (an approximation which should be valid at high densities) and normalize 
the off diagonal matrix elements as follows: 



The replicated free energy per particle, (j){m,T) = ^{m,T)/m, can be expanded in series: 

+ c+ log /3 r+ + c_ log /3 r_ + 

2m 



+ 



1 (m- 1) ^ /TrCP\ 



N 2m p 



where the p-th order term depends as usual on the p-points correlation function 

{TrCn = E E / d'z, . . . d% f^g'-Hzi • • • z,)Cl\%{z, - z,) ■ ■ ■ 

ei...epG{+,-} ^^l■■■^^P 

■ ■■CX-T.Mv-^ - ^v)CXk^^v - ^i), (23) 

We use a 'chain' approximation in the computations of traces, where terms with two equal 
indices are neglected , and the so called superposition approximation for the p-points cor- 
relation functions g^'''\zi...Zp) = g{zi — Z2) ■ ■ ■ g{zp — Zi). With these hypotheses we arrive 
at: 

(TrCO = / d'^z, ...d%ff E E a^'^'i^i - ^2)c;\%{z, -z,)--- 

fil...fip ei...ep 

■ ■■g'^-''^{zp-i - Zp)C;ii\X{zp^i - Zp)g^^^^{zp - z,)CX^^^{zp - z,). (24) 
The convolutions are computed in Fourier space, introducing the tensor: 
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D%{k) = I dW' {r)C%{T)e^^^ , (25) 

which can be decomposed into its diagonal (longitudinal) and traceless (transversal) parts 
with respect to the spatial (/x, v) indices: 

Di{k) = V a-\k) + - ^) h-\k). (26) 

The last step consists in the diagonalization of D in the space of the particles types (e, e'). 
For each k, there are four distinct eigenvalues, the two 'longitudinal' ones, corresponding to 
that of the 2x2 matrix 

Df\k) = a'^'ik) + ^6">), (27) 
and the two 'transverse' eigenvalues of the matrix 

Df'{k) = a'''{k)-h'''{k). (28) 

(Jj 

The eigenvalues are: 

All = I + D^- + ^(D++-Dp)2 + 4(D+-)2) 

i + Dp - ^(I)++-Dp)2 + 4(D+-)2 

l^^ = \ {d1+ + Dr ^ [Dt' - DTf + A{Dt-Y) (29) 

Using these approximations, the expression of the binary mixture free energy per particle is: 

, / X (i(m — 1) , , , dim — 1), , ,^ . , 
0(m, /5) = log(m) - log(27r) + (c+ log(/3r+) + c_ log(/?r_)) + 

+ ^^^^- / d''k{Ls{\\\{k)) + Lsiii\\ik)) + id-l) [L3(Ax(^)) + L3(M^))]} + 
- / d'rpj:9-'{r)Y:{C-:{r)f-^logZu,Wm), (30) 

where the function L3 is log(l — x) + a; + x'^/2. 

Let us notice that the condition for identifying the Kauzmann temperature, ^|^|m=i = 0, 
reads in our harmonic approximation: 

Su, = ^ log(27re) - ^(Tr log (/3M)) (31) 

Siig is the entropy of the liquid at the effective temperature T^jf, which is equal to T for 
m = 1. The right hand side of this equation is nothing but the entropy Sgoi of an harmonic 
solid with a matrix of second derivatives given by M. Thus, we find: 
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= m 



dm 



Sliq Ssol (^2) 



m=l 



If Siig < Ssol, the system is in the glassy phase (T < T^), while in the other case Sug > Ssoi, 
the temperature is greater than Tk (and of course less than To if the spectrum of M is 



positive). The complexity is then E = Sug — Ssoi, as expected on general grounds |T5|. 

Formula (|30D allows to compute the free energy $(m, T) = m(/)(m, T) which is the main 
quantity needed to investigate the thermodynamics of the low temperature glass phase, 
using (^). It should be emphasized that within the approximations we used here, the only 
properties of the liquid phase which are needed to get $ are the pair correlation g{r) and 
the free-energy. Beside usual thermodynamic quantities (energy, entropy, heat capacity...), 
we are interested in the two new parameters describing the glassy phase: 

• The square cage radii A^, defined as = |((a;f) — (xj)^) for type e particles. This 
square cage radii are obtained by differentiating the free energy with respect to coupling 
terms and by sending couplings to zero in the end: 

A 2 djPF) , 

d[m — 1)N^ c/(l/aj 

The square cage radii are nearly linear in temperature in the whole glassy phase, which 
is natural since non harmonic effects have been neglected. 

• The effective temperature T^jj = T/m oi the molecular liquid. This temperature 
varies very little and it remains close to the Kauzmann temperature when T spans the 
whole low temperature phase, confirming the validity of our description of the glass 
by means of a system of molecules remaining in the liquid phase. It is worth to stress 
that the linear behaviour of the parameter m as a function of T is a feature shared by 
every IRSB system to our knowledge. 

The harmonic expansion makes sense only if M has no negative eigenvalues, which is 
natural since it is intimately related to the vibration modes of the glass. Notice that here 
we cannot describe activated processes, and therefore we cannot see the tail of negative 
eigenvalues (with number decreasing as exp(— C/T) at low temperatures), which is always 
present It is known however that the fraction of negative eigenvalues of M becomes 



negligible below the dynamical transition temperature ||2^ . So our harmonic expansion 



makes sense if the effective temperature Tg/j is less than T^. 

B. Small cage expansion 

It is possible to introduce a slightly different way to compute the molecular liquid free-energy, 
in order to take into account: 

• Non-harmonic terms. 

• Corrections to the quenched approximation. 
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Starting from the expansion of the potential in powers of the relative variables m, if one 
expands also the exponential of the corrective term, one obtains an expansion of as a 
power series in and a_. This is the generalization to mixtures of the small cage expansion 
scheme utilized in the pure case . This expansion is not equivalent to computing pertur- 
batively quartic and higher order corrections to the Gaussian approximation represented by 
the harmonic re-summation. Indeed, in this a± expansion we are using a truncated version 
of the series in (^). On the other hand this direct expansion allows to take into account 
the annealed fluctuations of the matrix M which were neglected in the harmonic approxi- 
mation. Therefore these two types of approximations are complementary. In this paper we 
consider the harmonic re-summation and the small cage approximation as distinct schemes 
of approximation and we compare results obtained independently in both them. However, 
it is clear that a better approximation of the replicated free-energy could be obtained by 
adding corrections from the small cage approximation, treated in some systematic way, to 
the harmonic re-summation. A first attempt in this direction will be found in a following 



work 22 



The leading term of (|T6[) in the a+, a_ ^ oo limit is: 



/2 vr a+ 
m 



■dN+ (m-1) 



'2 TT a. 



m 



-dN- (m-1) 



(34) 



Accordingly, the zero-order free-energy is: 



/30^°^(a;+, a_, m, 13) = \- log = (Iq + Qq (c+ logQ;+ + c_ log 



niN 



a. 



(35) 



with 



do = — log — log m log Zug{P m) 

2m m 2m m JM 



d (m — 1] 



The first order term is: 



(36) 



1 /• / ^ \ / m iv 



vt=l 



m N 



\a=l i=l 



N 



nk'^(E<) 



i=l 



a=l 



exp 



4aj 



a,b jS+ 



a , b\2 



u: - - 



4Q!_ 



a,b ia- 



i<j 



i3 m d 

1 - f E E E« - - - z,) 

i<j a=l 



m d 



= zi^^ 1 - ^ ( E E E« - - ^%)v,Az^ - z,) 

from which we get the first-order contribution to the free-energy: 



(37) 
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f3(j)^^\a+, a_, m, f3) = c+ af a+ + c_ ai a_, 

where we define the first order coefficients as: 
, dim — 1 



(38) 



2m2 
d im — 1] 



d'rpg++{r)J2V++{r) + c 



2m2 



d r p g~ 



MM 



r + c. 



d'^'r p 
d'^r p g~ 



En 



MM 



MM 



(39) 



Up to first order, the harmonic re-summation and the small cage expansion give the same 
results. Differences appear at the second order level, which is presented in appendix. In 
fact, the second order term in the harmonic re-summation is: 



4m 



^-^ [d'rpY.9-'{r)Y.{Cli 

•> ee' M^^ 



(40) 



while the second order term in the small cage approximation adds two new kinds of term 
(see the Appendix): 

• Those involving fourth derivatives of potential, which are anharmonic corrections, are 
proportional to (m — 1)^, unlike any other term up to second order. This means 
that they are less important near Tx where m ~ 1, and more important at very low 
temperatures. 

• Those expressing the ffuctuations of the diagonal terms of M. These are corrections 
to the 'quenched' approximation. 

The free-energy per particle, up to second order, is then: 



/30(a+, a-,m, (3) = do -\ — - (c+ log + c_ log a_) + c+ ata+ + c_ a, a_ + 

m 



(41) 



where the coefficients a" are given in the Appendix. 

The free energy should be studied in the zero coupling limit, that is Q;+,a_ oo. 
This can not be done directly with a powers series of a+, truncated at a finite order. 
Therefore one must first take the Legendre transform of 0, as previously discussed, getting 
the thermodynamic potential G as an expansion in powers of different cage sizes A^, defined 
by means of (p3D. Within this formulation, the free energy in the vanishing coupling limit 
is obtained by looking for possible minima of G with respect to A'^,A~ . 
The Lagrange transformed free energy is, at first order: 



/5G(A+,A-,m,/?) =70 + 



d (1 — m) 



2m 

+ c+ 7^ A+ + c_ -f^A 
d{l — m 



c+ log(A+) + c_ log(A" 



7o = cto + 



m 



7i 



■-1 5 



7i 



■-1 5 



(42) 
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and the saddle points equations read: 

dG rf(l-m) 1 



dG ^ rf(l-m) 1 1 

- ^ A - 



dA~ m 'Ji 13 



- / d'r p gT-{r) \^V'-{r) + J d'r p g^;{r) ^AV+-(r). (43) 
The first order free energy in the vanishing coupling limit is correspondingly given by 

^G{A^\ A-*, m, /3) = ^il^ log - A log(^) _ -\ log ZU^m) + 

Im \m / 2m m IS 

- ^^^^^ (c+ log(/3 r+) + c_ log(/3 r_)) . (44) 

This expression for G looks quite reasonable. First of all one may note that in the m — >■ 
1 limit it reproduces the 'hquid' free energy density jSf = — log Ziig{P)/N, as it should. 
Moreover, in the limit in which the "+" and "— " particles are no more distinguishable, this 
expression for G coincides with the one found in the pure case (More precisely, the 
two generalized free energy would coincide in this limit if the liquid free energies at inverse 
temperature (3 m were the same, which would be true if one would forget about the mixture 
entropy contribution oc c+ logc+ + c_ logc_). 

The computation of the second order terms can be carried out in a very similar way. 
One gets (see the Appendix): 



PG{A+, A-,m, /3) = 70 + 73 (c+ log(A+) + c_ log(A-)) + 

+ c+ 7+ + c_ 7r A- + 

+ c+ 72++ (A^? + c- {A- f + c+c_7rA+ A-, (45) 

with 

In evaluating the formulae of the appendix one needs to know the three particles cor- 
relation function. This correlation function can be computed starting from a generalized 
HNC expansion Here we follow the simpler route of evaluating the three point function 
using the superposition principle, i.e. gs^x, y, z) = g{x — z)g{x — y)g{y — z). When looking 
for the minimum dGjdA^ = 0,dG/dA~ = 0, one faces the problem that the second order 
corrections are very important (this happens also in the pure case). In this case the solution 
can be found only through a perturbation around the first order solution. In this way one 
gets 

A+* = A+* + 6A+* 
A~* = A:[* + 6A^* 
G{A+*, A-*, m,(3) = Gi + 6 G2 (47) 
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and, by writing m = mi + Sm2, the stationarity condition reads: 



dm \ dm 



-2 = -i^K) (l^K) . (48) 



Therefore, one looks for the value ml which maximizes Gi, which is nothing but the first 
order free energy, and then one computes the second order corrections at m = mj. The 
result is: 



' ill)' ^ ill fit' 

G, = cJ^] ir + cJ^] 72- + c+c__^72+- (49) 
V7i / \li J li li 

Finally, the second order correction to mi is obtained by (^S]). 



IV. HNC FOR BINARY MIXTURES 



In evaluating the liquid free energy fug and the g^^' at the effective temperature T/m we 
use the so-called Hypernetted Chain approximation (HNC), a simple closure approximation 
that consists in neglecting the 'bridge' diagrams in the Mayer expansion p3H26|l. For homo- 
geneous fluids, apart from the constant Smisc = — c+log(c+) — c_ log(c_), the free energy of 
the liquid in the HNC approximation is given by: 

l/3F[{/^'(r)}] = logp - 1 + f / d'rY^c^c,, {^'(r) [log g^^' (r) + f3V^^' (r) - l] + l} + 



2pJ {2ttY 
where 

/i"'(r) 

and V is the determinant of the matrix 



(p h'''{k)f 



r ir) - 1, 



l + pc+h++{k) pc+h+-{k) 
pc_h+'{k) l + pc_h' {k) 



(50) 



(51) 



(52) 



The closed set of HNC equations for the two point correlations can be derived as a station- 
arity condition of the functional F with respect to these correlation: they are solved using 
the same numerical technique utilized in the pure case p[. 
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HNC is expected to be a good starting point for our study since both / and the mean 
values of quantities that involve only two particles correlation functions seem to be evaluated 
with an error smaller than 10% in the temperature region we are interested in, as we verified 
by comparing the analytical estimations with simulation results. The terms involving the 
three point correlation functions, when evaluated by the superposition approximation and 
the HNC pair correlations, are reproduced with errors which seem to be smaller than 30%. 

V. RESULTS AND DISCUSSION 

Before discussing the analytical and numerical results on the soft sphere binary mixture, 
let us pay attention, for a while, to the soft sphere model (P) with the particular value 
R = 1 (i.e., the pure case). This allows to compare thermodynamic quantities obtained 
within the small cage expansion up to second order (evaluating the three point function by 
the superposition approximation) with those computed at the same order in the replicated 
HNC re-summation scheme 

3.6 p 
3.4 - 
3.2 

3 - 

M- 2.8 - 

2.6 - 
2.4 - 
2.2 - 

^0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 

T 

FIG. 1. The free energy of the pure soft sphere model versus temperature. The three curves 
are the results obtained from the small cage expansion at the first order (dotted line) and at the 
second order (dashed line), and those from the HNC re-summation scheme (continuous line). 
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FIG. 2. (3 m (a) and [3 A (b) of the pure soft spheres model versus temperature. The three curves 
are the results obtained from the small cage expansion at the first order (dotted line) and at the 
second order (dashed line), and those from the HNC re-summation scheme (continuous line). 



In [Figg. l-2a,b] we show the free energy, the effective inverse temperature jSm and 
the cage radius A for the pure soft sphere model both at the first order, that gives the 
same results in the two cases, and at the second one. As already outlined, when starting 
from the generalized HNC expression, the second order coefficient 72 is obtained without 
further approximations than the one related to the use of HNC. On the other hand, these 
results confirm that evaluating the three point correlation function which appears in 72 
by the superposition approximation is a rather good approximation. In particular we get 
very similar values for the thermodynamic transition point, — 1.53 from the HNC 
resummation scheme and Tx — 1.49 when using small cage expansion, i.e. an error less 
than 3%. 

Now we come back to the soft spheres binary mixture with the interaction parameters 
described in (P), taking in particular the value R = 1.2 of the ratio between the effective 
diameters in order to obtain analytical results comparable to the numerical ones. We consider 
both the small cage expansion to second order and the harmonic re-summation, finding 
results in very good agreement as is shown in [Fig. 3], where the glassy phase free energy 
computed in the two different schemes of approximation is plotted as a function of T ( for 
simplicity we take in the following p = 1). 
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FIG. 3. Free energy of the soft sphere mixture vs temperature. The continuous hne is the result 
of the harmonic resummation scheme and the dashed hne is the result of the small cage expansion 
to second order. 

The evaluations of the thermodynamic critical temperature obtained by the two analytic 
methods nearly coincide: we get Tk = Tk'^ ~ 1.65, which is in agreement with the 
numerical estimates that we are going to discuss. For the sake of comparison, let us remember 
that the Mode Coupling critical value for this model is F^) ~ 1.45. Let us note that the 
ratio Td/Tk is usually found to be between 1.2 and 1.6. 

We stress that the parameter m and cages size, and A", plotted in [Figg. 4a,b] are 
nearly linear with temperature. This means, in particular, that the effective temperature 
T /mis always close to Tj^, so in our theoretical computation we need only the mean values 
of observables in the liquid phase, at temperatures where the HNC approximation still works 
quite well. 

One can also observe that the specific heat (see [Fig. 5b]) shows an evident 'jump' at T^, 
remaining close to the crystal-like value, 3/2 (we have not included the kinetic energy), in 
the whole glassy phase. The qualitative behavior of thermodynamic quantities, apart from 
the presence of the two distinct radii, is very similar to that observed in the pure case |]l],|^ 
and it corresponds to a second order transition from the thermodynamic point of view. 
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FIG. 4. In (a) we plot f3m vs temperature, from the harmonic re-summation scheme (continuous 
line) and from the low temperature expansion to second order (dashed line). In (b) we present PA~^ 
(continuous line) and PA~ (dashed line) computed in the low temperature expansion to second 
order. Note that, quite reasonably, the smallest cage radius corresponds to particles with the 
largest effective diameter. 




FIG. 5. The energy (a) and the specific heat (b) of the soft sphere mixture versus temperature, 
both in the liquid and in the glassy phase, from the harmonic re-summation scheme. 

The harmonic re-summation scheme suggests an intriguing approach for evaluating the 
thermodynamic critical temperature by simulations, starting from (|3^). Here the liquid 
entropy can be obtained for instance by numerically integrating the energy 

Su,{P) = P {Eu,{(3) - Fu,m = + P Eu,{(3) - f d(3'Eu,m (53) 

where 5*;°^ is the entropy of the perfect gas in the /5 — > limit, i.e in the binary mixture case 

Sl„ = iV (1 - log p - c+ log c+ - c_ log c_) . (54) 
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Moreover, one can think of directly numerically evaluating the 'harmonic solid' entropy 



= ^(1 + log(27r)) - ^ (Tr log(/3A^)) , (55) 

by diagonalizing the 'instantaneous' Hessian and by averaging over different configurations. 
The knowledge of Suq and Ssoi allows to obtain a numerical estimate of Tk as the temperature 
where the two entropies become equal, and to measure the complexity 

j:i(3) = ^[SH,i/3)-S,oim- (56) 
When attempting to obtain such evaluations, we face two kind of problems: 

• The well known hard task of thermalizing glass-forming liquids at low temperatures. 
Here we choose to perform a simulated annealing run of a quite large system, using data 
on the liquid energy down to the temperature where the equilibrium was still reachable 
in a reasonable CPU time (F ~ 1.5). Then we extrapolate the liquid entropy behavior 
at lower temperatures by fitting data in the interval F G [1, 1.5] with the power law 

SH,{T)=aT-^/' + b. (57) 



In fact, it has been shown ||30| that the potential energy of simple liquids at high 



densities and low temperature must follow this law, and we find that our numerical 
data are in a very good agreement with it. 

The correct evaluation of the solid entropy, which is a subtle task. Beyond the mean 
field approximation there always exists a non zero number of negative eigenvalues, 
which decreases as exp(— C/T) at low temperatures and is expected to be negligible 
below the Mode Coupling temperature. An estimate of the error on Sgoi can be found 
by doing the following two measurements, a) One includes in the computation of 
Trlog(/3A^) only the Mpos positive eigenvalues, b) One includes all eigenvalues, but 
one takes the absolute values of the negative ones: 



qia) 
^sol 

N 



S 



(b) 
sol 



N 



d 
2 

d 
2 



,27r 



■Afpos 



:i+iog(-T))-(Ti-Eiog^^) 



(58) 
(59) 



The percentage of non positive eigenvalues that we find by diagonalizing the 'instan- 
taneous' Hessian is still about 20% at F ~ 1, it decreases to less than 10% at F ~ 1.2 
and in the region definitely below T^, i.e. above F ~ 1.5, it is ~ 4%. On the other 
hand, nearly all the negative eigenvalues are less than one in absolute value. There- 
fore, particularly at temperatures T ~ To, we find a sizable difference between S*^"] 
and Sf^i (we disregard in both cases the very few |A| < 10~^), as is shown in [Fig. 6]. 
One should note that Sf^i seems to display the most regular behavior as a function of 
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r. The presence of negative eigenvalues is possibly related also to the fact that when 
diagonalizing the 'instantaneous' Hessian the system can be far from the 'center' of the 
minimum, in positions where higher order corrections to a harmonic approximation of 
the energy landscape are important. 

A more extensive study should be performed in order to better understand these sub- 
tleties of the computation of the solid entropy. However we would like to mention here 
a third way for evaluating numerically the solid entropy. Starting from an equilibrium 
configuration at a given F value, we performed a Monte Carlo run at T = 0, allowing 
only quite small displacements to each particle. The percentage of non positive eigen- 
values becomes very rapidly < 2% in the whole temperature range considered and 
correspondingly the two different ways of evaluating the 'solid entropy' give compat- 
ible results. The obtained S'^q] is near to the one evaluated from the 'instantaneous' 
Hessian by using also the absolute values of negative eigenvalues in the region T ~ To 
but it decreases slightly faster when lowering the temperature (see [Fig. 6]). 

More precisely, we performed a simulated annealing run of a system of iV = 258 particles, 
in a cubic box with periodic boundary conditions, starting from F = 0.05 and performing up 
to 2^^ MC steps at each AF = 0.05, the maximum shift 6max permitted to each particle in 
one step being chosen such as to get an acceptance rate ~ 0.5. The energy and its fluctuation 
were measured in the last half of the run at a given F- value. 

Just for decreasing the error on the evaluation of Suq, we fit the very high temperature 
data on the energy, up to Fq = 0.2, by using T^EiT) = aF^+6F+c, obtaining correspondingly 
F(Fo) = 4aFQ/3-|-26FQ + 4cFo that turns out to be perfectly compatible with the HNC value 
(i.e. we are still in the region where no differences are observable between numerical data and 
the HNC approximation). The integration is subsequently performed by interpolating with 
a standard numerical subroutine the simulation data in order to get a result independent on 
the integration interval. 

In order to evaluate S^^j and SfJi we considered 16 different configurations in the last 
half of the run at each F-value, while S'fJi was measured from the configurations obtained 
by these ones with 5000 MC steps at T = (starting from S^ax = 0.1 and decreasing it up 
to 0.02 during the run). One should note that at a given F-value the obtained evaluations 
of Sg'^Ji seem to depend weakly on the starting equilibrium configuration (i.e. fiuctuations 
are very small). Moreover we get perfectly compatible results both halving and doubling 
the number of MC steps (in the last case we find practically only positive eigenvalues) . 
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FIG. 6. The entropy of the hquid (+) and the different evaluations (see text) of the amorphous 
soUd entropy, s^"] (x), sf^^ (*) and s^^^ (□), as functions of F, obtained in the numerical simulation. 
In the liquid entropy case the line is the best fit to the power law suq = a0^/^+h, otherwise lines are 
only interpolations between neighboring points. The liquid and solid entropies seem to cross around 
Tk ~ 1.75, which is the corresponding estimation of the thermodynamic liquid-glass transition. 

We plot in [Fig. 6] both ^^^(r) and the obtained evaluations of Sgoiiy) by the different 
ways considered, s^^], s^g}i and s^^}i are very close to each other when approaching the 
thermodynamic liquid- glass transition, giving similar estimates of Vk ~ 1.75. 

The study of the system coupled to a reference configuration Xref-, which is an equilibrium 
configuration of the system itself at the considered temperature, allows to measure the 
complexity by an alternative route. One considers 

(5n = pHo + e{x- XrefT, (60) 

where 

(x-X,e/)^^EEC^r-<e// (61) 
i=l /x=l 

is the squared distance between the configurations (note that the coupling breaks the 
rotation-translation-permutation invariance) . Therefore 

/5/(e, (3) = PfM + r de'iix - Xreff),, (62) 
Jo 

where in the region T ~ T^, 

(3fo{(3) = lim^ /(e, (3) ~ (3e{(3) - E(/3) (63) 



22 



On the other hand, one has 



(64) 



Jim /?/(6,/?) = PfM = Pe{P) + ^ (log(^) - 
This means that one can obtain an evaluation of the configurational entropy as 

m ^ 4 + f de'iix - Xrefr)e' - ^ (log(^) - l) , (65) 

in the large e limit, taking into account as usual the perfect gas binary mixture entropy. 

Here wc considered a large system of = 2000 particles in a cubic box with periodic 
boundary conditions and we put a cut-off on the potentials, i.e. V"'' (r) = V" {Rmax) for 
r > Rmax, choosing Rmax — 1-7 that means a practically negligible V'^'^' {Rmax) ~ 10~^. The 
algorithm is then implemented in such a way that for each particle the map of the ones 
which are at distance lower than Rmax + '^^max is recorded and updated during the run. 

We performed up to M = 2^1 MC steps at each considered V = 1.4, 1.6, 1.8, 2.0. At 
the end of the run, the configuration was copied in the reference one and subsequently a 
run oiM/lB MC steps was performed on the coupled system for different e values, e = 1, 2, 
4, 8, ... up to very large e ~ 10^, measuring the squared distance. We note that perfectly 
compatible results were obtained for A/" = 2^^. The integrals were evaluated by interpolating 
with a standard numerical subroutine between the simulation data in order to obtain results 
independent on the integration interval. 
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FIG. 7. In (a) we plot {{x — Xref)^)e as a function of e at the considered values F = 1.4 (+), 
1.6 (x), 1.8 (*) and 2.0 (□). Here we show also 3/(2e) (the dashed line) in order to make ev- 
ident the reaching of the asymptotic behavior. In (b) we present the evaluations (see text) of 
S ~ s^jg -|- /q de'{{x — Xref)'^)e' — 3 (log(€/ (27r)) — 1) /2 in the large e limit (the different curves, from 
top to bottom, correspond to F = 1.4, 1.6, 1.8 and 2.0 respectively). The complexity turns out 
to be compatible with zero at Tk = 1.6 which is therefore the evaluation of the thermodynamic 
liquid-glass transition temperature. 
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Wc plot in [Fig. 7a] both the data on {{x — XrefY)e at different V as function of e and 
3/(2e). The asymptotic behavior seems to be reached around e = 2000, though also at 
larger e there are very weak deviations from it. When looking at the difference between 
the corresponding integrals and (i(log(e/27r) — 1) — s°j^ in the large e limit (see [Fig. 7b]), 
one finds that the complexity is compatible with zero at F^- ~ 1.6, a value slightly lower 
than the previously obtained F^ ~ 1.75 (the analytical estimation being Tk — 1.65). On 
the other hand, the 'errors' on these estimations are difficult to evaluate but they might be 
quite large, a more extensive numerical analysis being necessary in order to improve these 
results. 
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FIG. 8. The complexity S(r) computed in the harmonic re-summation scheme (continuous line) 
and the different numerical evaluations, i.e. suq — s^^] (+), suq — s^^^ (x), suq — s^^-J (*). The □ 
correspond to the S-values obtained by studying the coupled system at F = 1.4 and 1.6. 

At last we plot in [Fig. 8] the different numerical estimations of the configurational 
entropy E as a function of the temperature and the behavior obtained analytically in the 
harmonic re-summation scheme. In spite of the 'uncertainties' in the measures of S and in the 
analytical approximations (first of all related to the use of the HNC closure for evaluating 
liquid quantities), the agreement between theory and simulation looks quite satisfactory. 
We leave for future work both more extensive numerical studies and the improvement of the 
analytical results, that should allow a more careful comparison. 

ACKNOWLEDGMENTS 

The work of MM has been supported in part by the National Science Foundation un- 
der grant No. PHY94-07194. BC would like to thank the Physics Department of Rome 
University 'La Sapienza' where this work was partially developed during her PhD. 



24 



APPENDIX: SECOND ORDER COEFFICIENTS 



The second order expression of replicated partition function in small cage approximation is: 



/9 



^(2) = ^o )i_ti(yyv'^ 



)E«.-«.>)K-«P) + 



2, 



+ 



2-16 



E E ^^^(^^ - ^) EK - ^^)« 



(66) 



where distinction between terms which involve sums over particles of a given kind is required. 
Taking the logarithm of partition function, one finds the second order contribution to (f) 



/30^^^(q;_|_, q;_, /9) = c+ a^"^ + c_ 02 + c+ c_ a J q;_|_ • (67) 
The second order coefficients also depend on the three points correlation functions 

1 / 



(ri,r2) = 



E 5{xi-Xj-ri)5{xi-Xk-r2)). 



(68) 



One has 



d'r p g^^T) E (r) + - d'r p g-^'ir) ^ ^^.(r) 
c\ f d'^nd^r,p'g++^n,r,)Y.V;^-^{n)V;^-^{r,)+ 



4 m'^ 

_ f3^ (m - 1) 
4 

+ 2 c+ c_ / d'^n p'^++-(ri, r2) E C^(^i) ^"(^2) + 



+ 



+ 



(m - 1) 
T m3 



c+ / c?'^rp^++(r)^i;++(r) l7/(r)+ 



+ 



On = 



^ / ci-^r p^+-(r)^y+-(r) l^Vl^) 



4 

(m - 1) 
T m3 



c?'rp^+-(r)^i;;-,(r) + 



/d'^rp5+-(r)^V+-(r)V+-(r). 



(69) 



and the expression of 02 is obtained by changing the '+' in ' in the coefficient a^^ . 
To obtain the Legendre transform one must solve the system of linear equations 
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d{l/a+) 2 



d{l-m) 



dil — m) , , 
- o c_, (70) 



2 

and substitute the solutions into: 

n^r A+ A- n\ n\ d (l — Ul) d (l — ITl) A~ 



getting: 



/?G'(A+, A", m, /?) = 70 + 73 (c+ log + c_ log A') + c+ 7^+ A+ + c_ 7f A" + 

+ c+ 72"-^ (^+)' + c- 72" (^")' + c+c_72''~^^ • A- (72) 

with 

7o = Co - ao(l + log m)/m 73 ^ -a^jm 

ji — ■ji = m Gi 

(73) 

7++ = m^a^+ + m3(a]^+)2/(2 ao) 7^~ = m^as ~ + m3(ar~)2/(2 ao) 
7+- = m^at~. 
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